#include "ngstat.h" 

/*

ngstat (c) 2012-2013 Emanuele Raineri
emanuele.raineri@gmail.com


This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.


this program implements the formulae contained in
Luca Ferretti et al: 
Neutrality tests from high-throughput resequencing...
and
Luca Ferretti, Emanuele Raineri, Sebastian RO:

Neutrality tests for sequences with missing data

**************************************************/


void compute_all(int* pos, int ws, int **a,float **c, int *n, int **r ) {
    int x,s;
    float fh,fh_var,w,pi;
    for(x=0; x<ws; x++){
        fprintf(stdout,"pos:%d, nind:%d ",pos[x],n[x]);
        fh = fhat( x,  a, c,  n[x], &fh_var );
        fprintf(stdout,"fhat:%g fhat var:%g\n", fh, fh_var );
    }
    w=theta_w(a, n , ws,  r , &s);
    pi = theta_pi( r,a,c, n, ws );
    float var_pi_w, d =0;
    var_pi_w = variance_tajima(ws,3,  n,  r,  w );
    fprintf( stdout , "S=%d theta_w=%g theta_pi=%g ", s, w , pi );
    fprintf( stdout , "theta_pi-theta_w=%g var(theta_pi-theta_w)=%g D=%g\n", pi-w, var_pi_w,d );
}


int main(int argc, char* argv[]){
    splash();
    if (argc!=2) {
        fprintf(stderr,"usage: ngstat <window_size> \n");
        exit(1);
    }
   
    /*
    int **a : position,individual: genotype
    int **m : pos, individual : number of minor alleles
    int **r : pos, individual : read depth
    int *pos : genomic coordinate of each position in the window
    float **c : pos, individual : avg chromosome parsed in the position (max 2).
    */

    int  k,status,i,j,lc,lr,**a, **m, max_nind, *nind, ncols, *pos, *parsed,  **r , ws;
    float var, **c;
    FILE* mpileup;
    char *line,*fields[MAXFIELDS];
    
    line =malloc(MAXLINE);
    ws       = atoi(argv[1]);;
    parsed = malloc(PARSEFIELDS*sizeof(int));
    
    for (i=0;i<MAXFIELDS;i++){
        fields[i]=malloc(MAXLINE);     
    }
   
    
    //rows correspond to positions, columns to individuals
    
    a   = malloc(ws*sizeof(int*));
    r   = malloc(ws*sizeof(int*)); 
    m   = malloc(ws*sizeof(int*)); 
    c   = malloc(ws*sizeof(float*));
    
    
    nind   = malloc(ws*sizeof(int)); //number of individuals sequenced
    pos    = malloc(ws*sizeof(int));


    for(k=0; k<ws; k++) nind[k]=0;
    
    lc=0;lr=0; //line counters 
    while(1){
        
        fgets(line,MAXLINE,stdin);
        
        if (feof(stdin)) {
            if (lc>0)  compute_all( pos, lc , a, c, nind , r);
            break;
        }
        
        ncols=split_line(line,fields);
        nind[lc]= (ncols-4)/2;
        pos[lc]=atoi(fields[1]);

        
        r[lc]=calloc(nind[lc],sizeof(int));
        m[lc]=calloc(nind[lc],sizeof(int));
        a[lc]=calloc(nind[lc],sizeof(int));
        c[lc]=calloc(nind[lc],sizeof(float));



        for(k=0;k<nind[lc];k++){
            
            //parse each individual and fill r[k],m[k],a[k],c[k]
            for(j=0;j<PARSEFIELDS;j++) parsed[j]=0;
            parse_pileup(toupper(fields[2][0]),fields[4+k],fields[4+k+nind[lc]],parsed);
            

            r[lc][k]=parsed[1]+parsed[2]+parsed[3]+parsed[4];
            m[lc][k]=r[lc][k]-parsed[0];
	    c[lc][k]=est_seq_chr(r[lc][k]);
            a[lc][k]=est_gtype(m[lc][k],r[lc][k]);							
        
        }
        
        lc++; lr++;
        
        if ( lc == ws ) {
            compute_all(pos,lc, a, c, nind , r );
            lc = 0;
            }
    }

    fprintf(stderr,"I'm done after reading %d lines\n",lr);
    exit(0);

}
